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1.  Objective 


The  objectives  of  this  work  are  to  incorporate  dislocation  transport  into  a  continuum  crystal 
plasticity  formulation,  implement  the  model  into  a  finite  element  code,  and  evaluate  the  model’s 
ability  to  capture  effects  of  dislocation  aggregation,  specifically,  the  grain  boundary 
strengthening  effects  and  stress  concentrations  at  grain  boundaries. 


2.  Approach 


The  intent  is  to  develop  an  enhancement  to  traditional  continuum  crystal  plasticity  models  that 
can  be  readily  incorporated  into  existing  finite  element  implementations  and  run  with  minimal 
additional  computational  overhead.  This  goal  is  facilitated  by  targeting  explicit  dynamic 
solutions  where,  because  of  the  relatively  small  time  steps,  it  is  often  possible  to  impose 
additional  constraints  through  an  operator-split  approach.  Details  of  the  crystal  plasticity  model 
and  the  time  integration  are  reported  elsewhere  (7),  and  only  the  modifications  related  to  the 
constraint  of  the  dislocation  flux  are  presented  herein. 

2.1  Flux  Constraint 

The  basic  model  is  that  dislocations  move  along  the  slip  direction  from  one  element  to  another. 
The  flux  of  dislocations  crossing  an  element  boundary  for  each  slip  system  is 

0  =  Pais  vasa  •  nf  ,  (1) 

where  p%is  is  the  dislocation  density,  va  is  the  dislocation  velocity,  sa  is  the  slip  direction,  and 
rif  is  the  outward  normal  to  the  given  element  face.  Continuity  is  enforced  by  requiring  that  the 
dislocation  flux  exiting  through  an  element  face  equals  the  flux  entering  the  adjacent  element 
through  the  common  face.  The  dot  product  in  equation  1  accounts  for  the  orientation  of  the  slip 
system  with  respect  to  the  element  face.  It  is  zero  if  the  slip  direction  is  parallel  to  the  face.  This 
provides  the  opportunity  for  sharp  jumps  in  slip  rates  across  parallel  slip  planes  while  enforcing 
continuity  along  slip  planes. 

The  dislocation  density,  velocity  and  the  Burgers  vector,  b,  are  related  to  the  continuum  slip  rate 
by  Oro wan’s  equation 

Ya  =  Pdis  va  b.  (2) 

For  a  shared  element  face,  and  assuming  that  the  Burgers  vector  is  constant  and  that  the  slip 
directions  are  closely  aligned  across  the  interface,  the  continuity  error  in  the  accumulated 
dislocation  flux  between  elements  can  be  approximated  as 
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[(.YU) neighbor  ~  (7°%]  ’  «/  =  £/■  (3) 

The  subscripts  e  and  neighbor  on  the  accumulate  slip,  ya,  denote,  respectively,  the  element  of 
interest  and  the  neighbor  sharing  the  face.  The  approach  taken  is  to  penalize  the  error.  It  is  not 
necessary  to  assume  that  the  Burgers  vectors  are  equal  or  that  the  slip  directions  in  the  adjacent 
elements  are  aligned,  but  these  simplifications  are  made  for  expedience  and  coding  clarity  in  this 
initial  implementation. 

With  some  adjustments  (i),  the  driving  force  for  the  penalty  method,  equation  3  is  summed  as 

N  faces 

Z|  g1  |  ^  Y~Sign(ya)[(ya)neighbor  ~  (7°%]  |s“  •  ft/|  =  Ea.  (4) 

'  f'  f= 1  ^ 

Lf  is  the  distance  between  element  centroids  associated  with  the  particular  element  face.  It  is 
intended  to  provide  a  larger  penalty  if  the  slip  difference  occurs  on  a  smaller  spatial  grid.  Ea  in 
equation  4  is  a  weighted  sum  of  unsigned  slip  gradients  on  a  slip  system  for  a  given  element.  If 
the  value  is  positive,  slip  is  deficient  in  the  element  and  additional  slip  is  promoted.  If  the  value 
is  negative,  further  slip  is  impeded.  It  is  significant  that  a  constant  gradient  results  in  face 
contributions  that  sum  to  zero.  Hence,  a  constant  gradient  is  not  suppressed.  The  excess 
dislocations  associated  with  the  gradient  are  assumed  to  be  uniformly  distributed. 

Close  examination  of  equation  4  reveals  that  the  constraint  resembles  the  micro  force  balance  in 
several  gradient  formulations  (2,  5).  Hence,  with  a  small  coding  change  using  the  square  of  Lf  in 
the  denominator,  the  implementation  approximates  a  simplified  version  of  established  gradient 
models. 

The  flow  strength  in  the  power  law  slip  rate  model,  ga,  is  adjusted,  based  on  the  nonlocal 
contribution,  to  increase  or  reduce  the  resistance  to  continued  deformation 

ga  =  do  ~  PPi  tanh  (b2Ea  —  j  ■  (5) 

Here,  p1  and  p2  are  dimensionless  parameters,  p  is  some  average  shear  modulus  for  the  crystal, 
and  b  is  the  Burgers  vector.  The  base  crystal  strength,  g$,  could  be  a  function  of  slip  to  capture 
strain  hardening,  but  it  is  assumed  constant  in  these  analyses  to  simplify  interpretation  of  the 
results.  The  hyperbolic  tangent  function  is  used  to  cap  the  influence  of  the  nonlocal  term.  The 
contribution  is  approximately  linear,  ~gp2b2Ea,  while  the  argument  is  small,  and  it  is  caped  at  a 
constant  value  of  ±pp1.  The  penalty  parameters  are  chosen  as  p1  =  5  x  10-4  and  p2  =  10000. 
In  all  of  the  nonlocal  continuum  simulations,  the  base  flow  strength  is  g$  =  33.33  MPa;  the 
density  is  p0  =2.7  g/cm3;  the  Burgers  vector  is  b  =  0.3  nm;  and  shear  modulus  is 
p  =  25,650  MPa.  A  low  power  law  rate  exponent,  m  =  0.005,  is  chosen  to  provide  nearly  rate- 
independent  behavior  while  still  smoothing  the  transition  from  elastic  to  plastic  response  at  the 
slip  system  flow  strength,  ga . 
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2.2  Model  Geometry 

An  idealized  two-dimensional  crystal  geometry  is  used  for  these  analyses.  The  crystal  consists  of 
three  slip  systems  set  in  an  equilateral  triangle  configuration.  This  allows  a  multitude  of  slip 
modes.  The  model  was  implemented  in  the  large-scale  parallel,  explicit  finite  element  code 
ALE3D  (4),  and  details  are  provided  elsewhere  (/).  The  crystal  plasticity  constitutive  model 
existed  previously  (5)  and  the  strengthening  terms  in  equation  5  due  to  the  nonlocal  effects  were 
a  straightforward  addition. 

Constant  strain  quadrilateral  elements  with  hourglass  control  (6)  are  used  for  all  of  the 
simulations.  Elements  in  which  all  four  faces  are  adjacent  to  an  element  of  the  same  initial 
orientation  are  interior  elements,  and  gradients  are  computed  directly  as  indicated  in  equation  7. 
Elements  with  fewer  than  four  faces  contacting  regions  of  the  same  orientation  are  either  on 
grain  boundaries  or  model  boundaries.  The  face  is  flagged  for  these  elements,  and  a  parameter  is 
checked  to  see  whether  it  is  treated  as  a  zero  flux  boundary  or  a  free  boundary  with  no  slip 
impedance.  For  the  non-interior  elements  with  zero-flux  boundaries,  ghost  elements  with 
opposite  slip  are  assumed  across  the  flagged  faces.  Grain  boundaries  and  surfaces  with  applied 
boundary  conditions  are  treated  in  this  manner.  For  free  surfaces  and  periodic  boundaries,  the 
ghost  element  across  the  flagged  face  is  set  with  the  same  slip  so  that  these  interfaces  do  not 
contribute  to  the  gradient. 


3.  Results 


The  effect  of  the  slip  continuity  constraint  (slip  gradient)  is  evaluated  on  two  configurations, 
each  at  multiple  length  scales.  All  simulations  are  two-dimensional.  The  first  configuration  is 
simple  shear  of  a  single  crystal  with  one  of  the  slip  planes  initially  aligned  orthogonal  to  the 
shear  direction.  This  creates  single  slip  conditions  for  a  straightforward  evaluation  of  the  model. 
The  second  configuration  is  a  polycrystal  constructed  from  regular  hexagons.  The  orientation  of 
the  crystal  lattice  for  each  grain  is  chosen  at  random. 

3.1  Single  Crystal  Simulations 

Single  crystal  calculations  were  run  at  four  size  scales  using  a  20  x  100  mesh  of  square  elements 
(figure  1).  Velocity  boundary  conditions  are  applied  to  the  upper  and  lower  surfaces  to  shear  the 
top  of  the  crystal  to  the  right.  Initial  velocities  of  all  interior  nodes  are  prescribed  consistent  with 
simple  shear  to  eliminate  ringing  as  the  explicit  dynamic  calculation  starts.  Periodic  boundary 
conditions  are  applied  to  the  lateral  surfaces  to  mimic  an  infinitely  wide  crystal.  The  heights  of 
the  single  crystals  simulated  were  50,  5,  1,  and  0.5  pm,  and  the  width  of  the  simulation  box  was 
20%  of  the  height  in  each  case.  Although  the  width  is  irrelevant  with  the  periodic  boundary 
conditions,  multiple  elements  are  used  across  the  width  to  demonstrate  that  the  boundary 
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conditions  are  applied  properly.  Slip  transmission  is  restrained  on  the  upper  and  lower 
boundaries,  and  slip  transmission  is  unimpeded  on  the  lateral,  periodic  boundaries. 


Velocity  boundary 
condition 

Periodic  boundary 
conditions 


<r 


Lattice 

orientation 


Figure  1 .  Initial  configuration  for  the  single  crystal  and  the  crystal  lattice 

orientation.  The  bottom  is  fixed  and  the  top  is  moved  to  the  right. 
Periodic  boundary  conditions  are  applied  coupling  the  left  and  right 
hand  sides. 


The  nonlocal  strength  contribution  to  equation  5  is  shown  in  figure  2  for  the  four  crystal  sizes 
and  at  two  shear  strains.  The  largest  crystal  in  on  the  left  and  the  smallest  is  on  the  right.  The  top 
row  shows  the  distribution  at  a  shear  strain  of  0.03  and  the  shear  strain  is  0.05  in  the  bottom  row. 
The  color  scales  are  different  for  each  crystal  size,  but  the  scales  at  each  crystal  size  are  the  same 
at  0.03  and  0.05  shear  strains  to  highlight  the  evolution. 
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Figure  2.  Distribution  of  the  nonlocal  contribution  to  the  strength  on  the  slip  system  aligned  vertically  in  the  crystal, 
for  the  crystal  thicknesses  indicated,  at  shear  strains  of  0.03  and  0.05. 

The  magnitude  of  the  nonlocal  strength  contribution  and  the  relative  depth  that  the  distribution 
penetrates  from  the  surfaces  increases  as  the  crystal  thickness  is  decreased.  The  sharpest 
gradients  are  expected  at  the  crystal  surfaces  where  the  slip  transmission  is  impeded.  For  the 
50-pm  crystal  the  strength  is  increased  only  within  a  few  elements  of  the  surface,  corresponding 
to  a  few  microns.  The  central  portion  of  the  crystal  sees  no  gradient  or  strengthening  effect,  even 
as  the  strain  increases  from  0.03  to  0.05.  The  boundary  layer  also  appears  to  penetrate  a  few 
microns  in  the  5-pm-thick  crystal  simulation.  However,  a  smaller  portion  of  the  crystal  is  nearly 
free  of  gradient  effects  for  this  smaller  crystal.  At  yet  smaller  crystal  thicknesses,  the  gradient 
effect  penetrates  the  full  crystal  thickness,  and  the  evolution  with  increasing  deformation  is 
evident.  The  strength  at  the  center  is  elevated  significantly  by  a  strain  of  0.03,  and  it  continues  to 
increase  with  further  deformation. 
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The  minimum  and  maximum  values  are  indicated  in  each  of  the  plots.  Except  for  the  smallest 
crystal,  the  difference  between  the  maximum  and  minimum  increases  as  the  crystal  size 
decreases,  and  the  difference  also  increases  with  increasing  strain.  The  increase  with  strain 
indicates  that  the  gradient  is  still  evolving  at  a  shear  strain  of  0.05.  The  trends  are  different  for 
the  0.5-pm  crystal  because  the  gradient  strengthening  is  beginning  to  run  up  against  the  cap  set 
by  the  hyperbolic  tangent  function  in  equation  5.  With  continued  deformation,  the  gradient 
strengthening  is  becoming  more  uniform,  albeit  at  a  higher  level. 

The  momentum  balance  in  the  vertical  direction  (y-direction)  dictates  that  the  y-gradient  in  the  re¬ 
direction  stress  component  is  balanced  by  the  horizontal  (v-direction)  gradient  in  the  shear  stress. 
Since  the  periodic  boundary  conditions  require  that  all  horizontal  gradients  are  zero,  the  stress  in 
the  y-direction  should  be  constant.  The  magnitude  of  the  y-direction  stress  is  not  determined  by 
the  momentum  equation,  just  that  it  is  constant.  The  calculations  show  a  constant  y-direction 
stress  to  five  significant  digits. 

More  important  for  current  purposes  is  the  momentum  balance  in  the  horizontal  direction.  The 
lack  of  stress  gradients  in  the  horizontal  direction  requires  the  shear  stress  to  be  constant  through 
the  thickness.  The  simulations  show  that  the  shear  stress  is  constant  to  five  significant  digits. 
While  the  v-direction  stress  must  be  constant  in  the  v-direction,  the  symmetry  conditions  and 
momentum  equations  provide  no  further  constraints  restricting  its  gradient  in  the  y-direction. 

The  normalized  slip  rates  corresponding  to  the  configurations  in  figure  2  are  shown  in  figure  3. 
Again,  the  largest  crystal  is  on  the  left  and  the  smallest  is  on  the  right.  The  top  row  contains 
results  at  a  shear  strain  of  0.03,  and  results  at  a  shear  strain  of  0.05  are  shown  in  the  bottom  row. 
The  plots  are  normalized  by  the  applied  shear  strain  rate  so  that  a  value  of  1 .0  would  indicate  a 
uniform  shear.  Since  the  shear  stress  and  the  reference  strength  are  both  constant,  the  slip  rate  is 
approximately  related  to  the  gradient  term  in  equation  5.  Second-order  factors  influencing  the 
slip  rate  include  the  change  in  lattice  orientation  and  non-zero  components  of  the  v-direction  and 
y-direction  stresses  that  modify  the  resolved  shear  stress. 

The  slip  rate  follows  the  same  trends  as  the  nonlocal  hardening  contribution.  The  variation  in  slip 
rate  is  greater  for  the  smaller  crystals;  and,  except  for  the  smallest  crystal,  the  differences  are 
greater  with  increased  deformation.  As  a  result  of  the  boundary  conditions,  slip  rates  are  low  at 
top  and  bottom  boundaries  compared  to  the  center  regions.  This  results  in  the  sigmoidal 
deformed  shapes.  More  severe  differences  in  slip  rate  result  in  greater  deviation  from  a  linear 
shear  deformation  profile.  The  slip  rates  for  the  0.5-pm  crystal  become  more  uniform  at  the 
higher  deformation  because  the  gradient  term  is  capped  by  the  hyperbolic  tangent  function.  The 
strength  is  more  uniform,  which  results  in  a  more  uniform  slip  rate. 
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Figure  3.  Distribution  of  the  normalized  slip  rate  for  single  crystals  of  the  indicated  thickness  and  at  shear  strains 
of  0.03  and  0.05.  The  slip  rates  are  normalized  by  the  applied  shear  rate. 

The  shear  stress-shear  strain  responses  for  the  various  crystal  thicknesses  are  plotted  in  figure  4. 
The  curves  are  identical  through  the  linear  elastic  regime,  and  all  yield  at  approximately  the  same 
stress,  approximately  34  MPa.  The  applied  shear  rate  is  50  times  the  reference  shear  rate;  and, 
accounting  for  the  strain  rate  sensitivity,  the  apparent  yield  strength  is  calculated  from  the  power 
law  rate  equation  to  be  33.99  MPa  rather  than  the  reference  shear  strength  of  33.33  MPa. 
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Shear  Strain 


Figure  4.  Shear  stress-shear  strain  response  predicted  for  the  four 
crystal  thicknesses. 

The  stress  shows  the  expected  trend  of  increasing  strength  at  the  smaller  crystal  sizes.  This  is 
directly  related  to  the  nonlocal  strength  in  figure  2.  The  nonlinear  dependence  on  specimen 
dimensions  is  due  to  dividing  by  the  element  size  squared  in  the  modification  to  equation  5.  The 
strain  hardening  rate  is  fairly  consistent  with  increasing  strain  for  the  three  larger  crystal  sizes 
but  not  for  the  smallest.  Stress  in  the  0.5-prn  crystal  peaks  as  the  hyperbolic  tangent  function 
places  a  cap  on  the  nonlocal  strength  contribution.  This  is  also  consistent  with  the  results 
presented  in  figures  2  and  3.  A  final  observation  from  figure  4  is  the  kink  that  is  most  evident  in 
the  larger  two  specimens  near  a  strain  of  0.02.  This  marks  the  transition  from  single  slip  at  lower 
strains  to  slip  on  two  slip  systems  at  larger  strains.  As  the  crystal  lattice  rotates  and  stresses  build 
in  the  x  and  y  directions,  the  crystals  are  able  to  accommodate  the  deformation  more  easily  with 
multiple  active  systems.  Due  to  the  angle  of  the  slip  plane,  slip  constraints  at  the  boundaries  are 
not  as  severe  for  the  second  slip  system,  so  the  strain  hardening  rate  is  reduced. 

The  effect  of  mesh  resolution  on  the  solution  is  investigated  by  rerunning  the  5-pm-thick 
simulation  using  twice  as  many  elements  in  each  direction.  The  results  from  the  40  x  200  mesh 
are  shown  along  side  of  the  20  x  100  mesh  results  in  figure  5.  Other  than  the  expected 
differences  in  smoothness  of  the  fields,  the  nonlocal  stress  and  slip  rate  distributions  do  not 
appear  to  be  influenced  significantly  by  halving  the  mesh  size. 
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Figure  5.  Comparison  of  nonlocal  stress  and  the  slip  rate  for  the  20  x  100  and  40  x  200  element  simulations  of 
the  5-pm-thick  single  crystal.  The  results  are  shown  at  a  0.05  shear  strain. 

For  a  more  quantitative  assessment,  the  applied  shear  stress  computed  at  the  crystal  boundaries  is 
2.3%  lower  for  the  finer  mesh.  Part  of  the  difference  may  be  due  to  the  lower  energy  solution 
expected  with  an  increased  number  of  degrees  of  freedom,  and  the  remainder  can  be  attributed  to 
discretization  error  associated  with  the  nonlocal  computations  and  the  operator  split  algorithm. 
The  time  step  was  also  a  factor  of  two  lower  in  the  fine  mesh  calculation  due  to  the  dependence 
of  the  Courant  stable  time  step  on  the  mesh  size.  The  smaller  time  step  improves  the  accuracy  of 
the  operator  split  integration. 

The  simulation  of  the  0.5-pm-thick  crystal  experienced  numerical  instabilities  when  the  strain 
increment  per  step  was  too  great.  This  is  thought  to  be  associated  with  the  operator  split  where 
the  strength  increase  from  the  slip  gradients  creates  a  driving  force  that  is  too  large  and  over¬ 
corrects  the  slip  rate.  Specifically,  with  a  strain  increment  of  1.54  x  10  7  per  time  step,  the  slip 
rate  for  the  0.5-pm  crystal  was  erratic  and  non-zero  only  in  scattered,  isolated  elements.  These 
isolated  regions  of  slip  occurred  briefly  and  died  out  quickly  as  deformation  proceeds,  and 
eventually  strain  was  incremented  in  the  entire  domain,  albeit  unevenly.  When  the  time  step  was 
reduced  by  a  factor  of  two,  such  that  the  strain  increment  per  step  was  7.43  x  10-8,  the 
calculation  was  well  behaved.  The  results  in  figures  2  through  4  were  run  with  a  strain  increment 
of  3.853  x  10  8  to  be  certain  that  the  time  step  was  small  enough  to  suppress  the  instability.  The 
calculations  for  the  larger  crystal  sizes  experienced  less  gradient  hardening,  and  they  were  run  at 
the  Courant  stable  time  step  without  any  additional  time  step  controls. 
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3.2  Polycrystal  Simulations 


An  idealized  polycrystal  was  created  by  filling  a  rectangular  region  with  regular  hexagonal 
grains  (figure  6a).  Simple  shear  boundary  conditions  were  applied  by  prescribing  velocities  to 
the  nodes  on  the  upper  and  lower  surfaces.  Periodic  boundary  conditions  were  applied  on  the  left 
and  right  surfaces.  The  orientation  of  the  triangular  crystal  lattice  in  each  of  the  grains  was 
random,  and  the  rotation  angle  from  the  reference  orientation  is  indicated  on  the  plot.  The  half 
grains  at  the  same  height  on  the  left  and  right  of  the  model  region  were  given  the  same 
orientation  to  facilitate  application  of  periodic  boundary  conditions.  As  with  the  single  crystal 
simulations,  the  initial  velocity  of  all  interior  nodes  was  prescribed  to  eliminate  ringing  from 
abrupt  imposition  of  boundary  conditions. 
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Figure  6.  Grain  structure  (a)  and  finite  element  mesh  and  (b)  for  the  polycrystal  simulations. 

The  default  inter-element  slip  rate  condition  for  all  elements  is  that  any  element  face  touching 
another  grain  will  have  restricted  slip.  This  is  imposed  by  assuming  that  a  ghost  element  across 
the  interface  has  equal  and  opposite  slip  in  equation  4.  The  restricted  slip  condition  is  enforced 
on  the  upper  and  lower  surfaces  and  on  grain  boundaries,  including  those  grain  boundaries  on  the 
periodic  surfaces.  The  half  crystals  on  the  periodic  surfaces  are  treated  differently;  the  element 
across  the  interfaces  is  assumed  to  have  the  same  slip.  This  is  not  a  truly  periodic  condition,  but  a 
data  structure  identifying  periodic  neighboring  elements  is  not  yet  available. 

Three  model  sizes  are  investigated:  160V3  /tm  x  300  /tm,  16V3~/rm  x  30  /tm,  and 
1.6V3~/um  x  3  /urn.  All  use  the  same  mesh  configuration,  scaled  to  give  the  appropriate 
dimensions.  Each  of  the  88  hexagons  was  discretized  by  an  identical  mesh  of  21 12  quadrilateral 
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elements  (figure  6b).  Nodes  are  shared  along  the  grain  boundaries,  so  the  deformation  is 
continuous  throughout. 


Slip  rates  normalized  by  the  applied  shear  rate  are  shown  in  figure  7  for  the  three  polycrystal 
sizes.  At  the  larger  crystal  size  the  strain  rate  localizes  into  well-defined  bands.  The  majority  of 
the  deformation  is  carried  by  two  horizontal  bands  with  some  scattered  activity  in  the  central 
region  of  the  model.  Blocks  of  grains  appear  to  remain  elastic  while  localized  shear  along  the 
horizontal  and  vertical  bands  accommodates  the  deformation  between  neighboring  blocks.  The 
bands  do  not  follow  the  grain  boundaries,  but  many  are  associated  with  grain  boundary  triple 
points. 


Figure  7.  Normalized  slip  rate  distribution  for  simple  shear  deformation  of  idealized  polycrystals  with  heights  of 
330,  33,  and  3.0  pm.  The  color  scale  is  the  same  for  all  three  plots. 

At  the  intermediate  polycrystal  size  the  nonlocal  slip  constraint  diffuses  the  deformation.  The 
slip  bands  are  still  fairly  well  defined,  but  the  peak  strain  rates  are  not  as  high  and  regions  of 
nearly  elastic  behavior  are  smaller  and  less  well  defined.  The  slip  rates  in  the  3.3-pm-thick 
polycrystal  are  considerably  more  diffuse  and  the  material  near  the  highly  constrained  top  and 
bottom  boundaries  has  the  lowest  strain  rates.  Grain  outlines  are  evident  as  the  slip  rate  tends  to 
be  high  or  low  at  the  grain  boundaries,  and  the  color  contrast  across  the  boundaries  accentuates 
them. 

The  nonlocal  strengthening  associated  with  the  gradients  is  shown  in  figure  8  for  all  three  slip 
systems  and  the  three  crystal  sizes.  The  color  scale  in  each  row  is  the  same  so  that  the  magnitude 
of  the  effect  of  the  slip  systems  can  be  compared.  The  scales  are  different  for  each  crystal  size  as 
the  strengthening  is  much  greater  in  the  smaller  model  region.  The  color  scale  for  the  330  mm 
polycrystal  is  set  to  a  relatively  low  value  of  0.5  MPa,  and  even  then,  the  gradient  contribution  is 
only  evident  at  the  grain  boundaries  or  near  the  most  highly  shear  regions  shown  in  figure  7.  The 
impact  on  the  solution  is  minor. 
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Figure  8.  Nonlocal  strength  contributions  on  the  three  slip  systems  for  polycrystal  model  sizes  of  330,  33,  and 
3.3  pm.  The  color  scales  are  consistent  within  each  row. 

The  nonlocal  strength  distribution  in  the  3 3 -pm  polycrystal  appears  in  the  grain  interiors  as  well 
as  at  the  grain  boundaries.  The  strong  interior  features  are  associated  with  stress  concentrations 
at  grain  boundary  triple  points,  and  most  correspond  to  elevated  slip  activity  in  figure  7.  Many  of 
the  grain  boundaries  show  strengthening  on  one  side  and  softening  on  the  other.  These 
correspond  to  increasing  slip  when  approaching  grain  boundaries  for  strengthening  and 
decreasing  slip  when  approaching  the  boundaries  for  softening. 
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The  nonlocal  strengthening  is  quite  prominent  in  the  3.3-pm  polycrystal.  As  with  the 
intermediate  size  polycrystal,  the  sign  of  the  gradient  effect  is  often  flipped  across  the  grain 
boundaries.  For  most  grains,  the  gradient  is  strongest  at  the  grain  boundaries  and  decays  toward 
the  grain  center.  However,  there  are  a  few  notable  grains  where  the  peak  values  are  in  the 
interiors.  These  correspond  to  locations  of  intersecting  slip  activity  in  figure  7.  The  hyperbolic 
tangent  function  causes  the  gradient  effect  to  saturate  at  a  level  between  12  and  13  MPa.  It  is  also 
notable  that  the  gradient  strengthening  occurs  on  only  two  of  the  three  slip  systems.  This  reflects 
the  lack  of  redundant  slip  for  the  idealized  crystal.  Only  two  slip  systems  are  active  at  any  time. 

The  shear  stress-shear  strain  response  for  the  three  polycrystal  sizes  is  presented  in  figure  9.  As 
with  the  single  crystals,  the  stress  is  higher  for  the  smaller  polycrystals.  The  nonlinearity  with 
length  scale  is  also  clearly  evident.  There  are,  however,  two  important  distinctions  from  the 
single  crystal  results.  The  first  is  that  the  initial  yield  point  varies  with  crystal  size,  whereas  it  did 
not  for  the  single  crystal  simulations.  This  is  thought  to  be  related  to  the  single  crystals  yielding 
throughout  simultaneously  whereas  the  polycrystal  yields  gradually  and  may  build  up  local  slip 
gradients  before  the  macroscopic  yield  is  evident. 


ShearStrain 


Figure  9.  Shear  stress  strain  response  for  three  different  size 
scales  of  idealized  polycrystals. 

The  other  notable  difference  is  that  the  curves  are  not  smooth.  This  could  result  from  a 
combination  of  the  evolution  of  the  crystal  lattice  orientation,  evolution  of  the  slip  gradients,  and 
wave  propagation  in  the  explicit  dynamic  calculation.  The  change  in  lattice  orientation  is  shown 
in  figure  10  for  the  largest  and  smallest  size  scales.  In  the  330-pm  polycrystal,  where  the  strain 
localization  is  more  pronounced,  lattice  reorientation  is  also  localized.  The  local  geometric 
softening  facilitates  the  shear.  In  contrast,  due  to  the  slip  continuity  and  gradient  constraints,  the 
lattice  reorientation  in  the  3.3-pm  polycrystal  is  smoothed  over  a  larger  region  relative  to  the 
grain  size,  and  the  lattice  within  a  grain  rotates  nearly  uniformly.  A  larger  portion  of  the 
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polycrystal  had  to  realign,  which  takes  longer  and  results  in  a  greater  load  excursion  before  it 
settles  into  a  nearly  steady  shear  mode. 


Figure  10.  Change  in  crystal  lattice  orientation,  in  degrees,  at  0.025  shear  strain  for  the  330-  and  3.3-pm-high 
polycrystals. 


4.  Physical  Size  Scale  Considerations 


The  motivation  for  imposing  dislocation  flux  constraints  is  that  dislocations  are  discrete  entities 
that  propagate  from  one  element  to  the  next  as  part  of  the  slip  process.  Another  aspect  of  the 
dislocation  discreteness  is  their  spacing,  which  is  typically  quantified  in  terms  of  the  dislocation 

C.  H  _ O 

density.  For  well-annealed  metals,  a  typical  dislocation  density  is  10  -  10  cm  ;  at  a  few 
percent  deformation,  this  increases  to  10s  -  109  cm-2,  and  for  a  very  heavily  deformed 
polycrystal,  the  dislocation  density  is  in  the  neighborhood  of  1011  cm  2  (7).  Table  1  lists  the 
element  areas  for  each  of  the  four  single  crystal  simulations  and  the  higher  values  of  dislocation 
density  for  well-annealed,  lightly  deformed  and  heavily  deformed  polycrystalline  metals.  From 
these  values,  an  average  number  of  dislocations  enclosed  by  an  element  is  calculated. 


Table  1.  Average  number  of  dislocations  per  element  expected  in  the  simulations  of  for  well-annealed,  lightly 
deformed  and  heavily  deformed  metals. 


50  jum  Crystal 
Elength  =  0.5pm 
Earea  =  0.25pm2 

5  jum  Crystal 
Elength  =  0.05pm 
Earea  =  0.0025pm2 

1  jum  Crystal 
Elength  =  0.01pm 

Earea  =  0.0001pm2 

0.5  jum  Crystal 
Elength  “  0.005jim 
Earea  =  0.000025jLim2 

107  cm  2  (10  1  pnT2) 

1/40 

1/4000 

1 0  5 

2.5  x  ltr6 

109  cm  2  (101  pm  2) 

2.5 

1/40 

1/1000 

1/4000 

1011  cm  2  (103  pm  2) 

250 

2.5 

1/10 

1/40 
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The  entries  in  table  1  that  are  less  than  1 .0  indicate  that  not  every  element  will  contain  a 
dislocation.  For  example,  1/4000  means  that  one  out  of  every  4000  elements  can  be  expected  to 
contain  a  dislocation.  An  implicit  assumption  in  continuum  crystal  plasticity  models  is  that  the 
dislocation  content  in  the  elements  is  sufficient  for  slip  to  be  smooth  and  continuous.  It  is  clear 
that  these  conditions  are  not  met  for  a  well-annealed  metal  using  any  of  the  single  crystal  meshes 
since  not  every  element  would  contain  even  one  dislocation.  Using  the  discretization  provided  by 
the  50-pm-thick  crystal  simulation,  a  sufficient  number  of  dislocations  would  be  represented 
within  each  element  when  the  crystal  is  heavily  deformed,  but  not  in  the  deformation  leading  up 
to  that  state.  Considering  that  dislocations  are  usually  not  uniformly  distributed,  element  sizes  of 
a  few  microns  may  be  necessary  to  assure  a  sufficient  number  of  dislocations  per  element  for  a 
proper  continuum  crystal  model  representation. 

A  further  concern  is  that  dislocations  are  usually  less  uniformly  distributed  in  heavily  deformed 
metals.  They  typically  organize  into  walls,  which  create  a  cell  structure.  The  cell  walls  have  a 
very  high  dislocation  density,  and  the  cell  interiors  have  a  low  dislocation  density.  Extensive 
analysis  of  heavily  deformed  nickel  by  Hughes  and  Hansen  ( 8 )  shows  that  these  cell  sizes  are 
greater  than  0.1  pm  after  cold  rolling  to  a  98%  reduction,  or  a  logarithmic  strain  in  excess  of  3.9. 
In  order  to  have  a  smeared  representation  of  such  a  cell  microstructure  within  each  element,  the 
element  size  would  have  to  be  on  the  order  of  1  pm.  If  the  elements  are  small  enough  to  resolve 
the  cell  structure,  there  are  discrete  lattice  orientation  jumps  across  the  cell  walls,  not  smooth 
gradients.  The  crystal  plasticity  models  will  have  to  be  enhanced  in  an  alternative  manner  to 
include  the  additional  deformation  mechanisms. 

4.1  Semi-discrete  Dislocation  Model 

In  an  attempt  to  push  continuum  finite  element  simulations  to  smaller  length  scales  where 
dislocations  are  sparse  within  the  elements,  a  semi-discrete  model  was  developed.  It  is  run  within 
a  standard  explicit-dynamic  finite  element  framework  that  is  described  in  the  first  year  DRI 
report  (9).  The  single  slip  constitutive  model  follows  the  traditional  formulation  with  three  major 
modifications:  (1)  only  elements  that  contain  dislocations  or  dislocation  sources  can  slip;  (2)  the 
slip  increment  is  quantized  in  terms  of  the  Burgers  vector  and  element  size;  and  (3)  elements 
designated  to  contain  dislocation  nucleation  sources  have  a  reduced  flow  strength.  In  addition  to 
these  slip  model  modifications,  the  code  tracks  dislocations  moving  from  one  element  to  another, 
and  it  also  tracks  the  total  number  of  dislocations  that  have  traversed  each  element.  Details  of  the 
model  and  complete  results  have  been  reported  elsewhere  (i). 

Results  from  a  two-dimension,  simple  shear  simulation,  50  pm  wide  and  100  pm  high,  are 
presented  in  figures  1 1  and  12.  The  mesh  was  400  x  800  elements,  for  an  element  size  of 
0.125  pm.  There  were  158  dislocation  nucleation  sites  distributed  randomly  in  the  mesh, 

/u  _ o 

providing  a  nucleation  site  density  of  approximately  3.2  x  10  cm  .  This  is  in  the  range  of  the 
dislocation  density  for  a  well-annealed  metal. 
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Figure  1 1 .  Number  of  dislocations  passing  through  each  element  for  the  discrete  dislocation  simulations.  The  center 
80  pm  is  omitted  from  the  100-pm  crystal  to  highlight  the  gradients  at  the  top  and  bottom  surfaces. 


Figure  12.  Number  of  dislocations  currently  within  each  element  for  the  discrete  dislocation  simulations.  The  center 
80  pm  is  omitted  from  the  100-pm  crystal  to  highlight  the  gradients  at  the  top  and  bottom  surfaces. 

The  number  of  dislocations  that  passed  through  each  element  at  a  shear  strain  of  0.01  is 
presented  in  figure  11.  Only  the  top  and  bottom  10  pm  of  the  100-pm  model  are  shown  since 
most  of  the  center  section  appears  as  lines  connecting  the  upper  and  lower  portions.  The  most 
obvious  feature  is  the  discrete  deformation.  Slip  occurs  only  along  slip  systems  containing  the 
nucleation  sites.  The  slip  traverses  the  crystal  vertically  along  lines  of  elements  that  contains  the 
slip  planes.  The  remaining  elements  are  elastic.  An  important  feature  of  figure  1 1  is  the  slip 
gradient.  Since  dislocations  cannot  pass  through  the  upper  and  lower  boundaries,  the  slip 
(number  of  dislocations  passed)  at  these  surfaces  is  zero.  The  greatest  number  of  dislocation  has 
passed  near  the  center  of  the  crystal. 

All  of  the  dislocations  must  lie  between  the  nucleation  sites  where  they  originate  and  the  crystal 
boundaries.  The  slip  distribution  is  directly  related  to  the  current  positions  of  the  dislocations. 
The  dislocation  positions  are  shown  in  figure  12.  As  with  figure  11,  only  the  upper  and  lower 
10  pm  are  shown  for  the  100-pm-thick  crystal.  Dislocations  of  opposite  sign  originating  at  the 
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dipoles  are  shown  by  the  red  and  blue.  Dislocations  of  one  sign  move  to  the  top  and  those  of  the 
other  sign  move  to  the  bottom.  The  dislocation  density  is  greatest  at  the  boundaries  and  tapers  off 
toward  the  center  of  the  crystal.  This  is  the  classic  picture  of  an  edge  dislocation  pile-up.  The 
dark  red  and  blue  elements  on  the  crystal  interior  indicate  the  location  of  the  dipole  nucleation 
sites  where  dislocations  can  accumulate  before  gliding  toward  the  boundaries. 

Each  element  can  contain  multiple  dislocations,  with  the  maximum  number  being  four  for  this 
simulation.  The  gradient  is  evident  in  the  sparseness  of  elements  containing  dislocations.  Near 
the  surface,  most  elements  along  the  slip  planes  contain  dislocations,  but  as  the  distance  from  the 
boundary  increases,  the  elements  containing  dislocations  become  increasingly  sparse. 

While  the  results  of  this  simple  semi-discrete  model  show  dislocation  pile-ups  and  other 
expected  features  not  captured  by  standard  continuum  models,  the  specific  approach  is  will  not 
produce  meaningful  results  at  smaller  size  scales.  The  overriding  issue  is  the  singular  stress  field 
from  the  dislocations  dominating  the  solution  as  the  mesh  is  refined  to  resolve  the  stress 
gradients.  The  stress  field  is  driven  by  the  quantization  of  slip,  which  provides  a  higher  stress 
magnitude  related  to  the  better  resolution  of  the  singularity  at  finer  spatial  resolutions.  The 
unresolved  singularity  dominates  and  pollutes  the  solution.  There  does  not  appear  to  be  a  range 
of  element  sizes  that  will  both  resolve  the  stress  field  and  not  suffer  from  the  effects  of  the 
singularity.  Perhaps  other  semi-discrete  approaches  could  be  successful;  this  one  was  not. 

Discrete  dislocation  dynamics  simulations  (10,  11)  explicitly  account  for  dislocations  and  their 
interactions  and  provide  one  means  for  incorporating  dislocation  microstructure  at  finer  spatial 
resolutions.  Finite  element  methods  have  been  coupled  with  the  discrete  dislocation  simulations 
though  several  approaches  (12-14).  These  types  of  formulations  should  be  employed  at  the  finer 
spatial  resolutions  in  multiscale  modeling  schemes.  There  could  be  a  transition  from  traditional 
continuum  crystal  plasticity  to  proper  discrete  representation  when  the  spatial  resolution  is  fine 
enough  (15). 


5.  Conclusions 


A  nonlocal  crystal  plasticity  model,  motivated  by  slip  continuity  between  neighboring  finite 
elements,  was  implemented  in  a  large-scale  parallel  finite  element  code.  This  is  the  first  large- 
scale  explicit  implementation,  and  the  first  solution  illustrating  gradient  effects  at  grain 
boundaries  in  a  large,  dynamic  polycrystal  simulation.  The  results  show  the  expected  trends  of 
decreasing  the  severity  of  gradients  and  increasing  strength  with  decreasing  physical  size.  In 
terms  of  the  original  objectives,  the  program  was  successful. 

It  is,  however,  a  qualified  success.  The  gradient  effects  are  only  significant  on  the  length  scale  of 
approximately  100  of  pm  or  less.  It  was  determined  that  spatial  discretization  for  such  problems 
can  be  on  the  order  of,  or  smaller  than,  the  scale  of  the  underlying  microstructure.  For  lightly 
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deformed  metals,  only  a  small  fraction  of  the  elements  would  contain  dislocations,  and  the 
continuum  crystal  plasticity  model  will  not  apply.  For  heavily  deformed  metals  with  a 
dislocation  cell  structure,  the  element  size  can  be  on  the  order  of  the  cell  structure.  The  discrete 
orientation  jumps  and  deformation  mechanisms  associated  with  the  dislocation  cells  are  not 
adequately  represented  by  traditional  crystal  plasticity  or  the  gradient  model. 

Hence,  while  the  model  can  be  run  to  obtain  results  at  small  size  scales,  the  solutions  will  not 
capture  the  additional  physical  mechanisms  associated  with  microstructure  at  the  enhanced 
resolution.  The  model  will  not  provide  significant  insight  into  deformation  at  the  sub-micron 
scale  since  the  model  does  not  adequately  represent  the  structure  or  mechanisms  at  the  sub¬ 
micron  scale. 
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7.  Transitions 


An  ARL  Technical  Report  (ARL-TR-6307)  with  a  more  complete  description  of  the  models  and 
results  is  in  press.  Results  of  the  model  limitations  were  discussed  at  the  Materials  in  Extreme 
Dynamic  Environments  Fall  Meeting. 
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